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Abstract 

o 

Q | We have calculated the first and second sound modes of a dilute interacting 

Bose gas in a spherical trap for temperatures (0.6 < T/T c < 1.2) and for 
systems with 10 4 to 10 s particles. The second sound modes (which exist only 
' below T c ) generally have a stronger temperature dependence than the first 

sound modes. The puzzling temperature variations of the sound modes near 
T c recently observed at JILA in systems with 10 3 particles match surprisingly 



o 
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well with those of the first and second sound modes of much larger systems. 



Since the discovery of Bose-Einstein condensation in atomic gases of alkali atoms |I] , there 



has been great interest in the broken gauge symmetry (i.e., the "phase") of the condensate. 
In the case of 4 He, its "phase" dynamics leads to the existence of second sound, which 

X" 

is essentially the out of phase pressure and temperature oscillations. In a series of sound 
experiments, Jin et.al at JILA [0 have observed a number of "puzzling" behaviors in the 
temperature dependence and the dissipation of the sound modes above 0.5T C . There are no 
explanations for these behaviors so far. Jin et.al. have speculated that the observed "m = 0" 
mode could be the "second sound". If this were true, it would be a demonstration of the 
broken gauge symmetry of the system. However, in the absence of a detailed calculation 
consistent with experiments, the identification of the second sound mode would be difficult. 

To help identify the nature of the sound modes, we have solved the linearized two-fluid 
hydrodynamic equations of an interacting dilute Bose gas in a spherical harmonic trap. 
We have in mind systems that are sufficiently large so that the hydrodynamic approach is 



accurate jp. It should be noted that the recent experiments at JILA |2| were performed 
on small systems with a few thousand atoms. While the hydrodynamic modes of a large 
system may be different from the sound modes of a small one, the study of the former is 
important in its own right. After all, the number of atoms in the Bose condensate has 
increased from 10 3 to 10 6 within six months after the initial discovery [|I]]. It would not be 
surprising if Bose condensates with 10 9 atoms were produced in the near future. On the 
other hand, the hydrodynamic modes of large systems are relevant for the sound modes 
of small ones, as the former must evolve smoothly into the latter as the number of atoms 
is decreased continuously. This suggests the possibility of identifying the nature of sound 
modes of a small system by studying their hydrodynamic counterparts in a large one. Indeed, 
comparing our results (for systems with 10 4 to 10 8 atoms) with the JILA observations 0, we 
find that the temperature variations of the observed sound modes show up in the analogous 
modes of the larger systems in a spherical trap. In particular, the "mysterious" behaviors 
of the JILA (m = 0) and (m = 2) modes in the range 0.5 < T/T co < 0.8 match closely 
with the behaviors of the second sound modes of the larger systems in same temperature 
range, while the frequency and temperature dependence of the observed (m = 0) mode 
above T c are identical to those of the first sound mode in the same temperature range. (The 
temperatures T co and T c are the transition temperature for the ideal Bose and the dilute 
interacting Bose gas respectively.) 

Our choice of spherical symmetry is to keep the calculations manageable. Moreover, as 
a first step, we shall ignore dissipation. While it is entirely feasible within our scheme to 
include dissipative effects, we feel that it is important (as in bulk 4 He) to first understand 
dissipationless hydrodynamics, so that one can clearly identify the dissipative effects later 
in a complete solution. 

Linearized Hydrodynamics : We begin with the two-fluid hydrodynamic equations 
of Bosons with mass M in an external potential <f)(r) |J, Mh = — Vg, gi = — 7iVj0 — VjlLj, 
s = — V • (sv n ), and v s = — -gV(/i + + Mv n • v s ), where n, g, LTy, s, /i are the number 
density, the momentum density, the stress tensor, the entropy density, and the chemical 
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potential respectively. Here, v n and v s = (h/M)V6 are the normal fluid and superfluid 
velocities respectively, where 9 is the phase of condensate. For a spherical harmonic trap 
with frequency lot, 4>{r) = \Moj\r 2 . In the presence of a condensate, n, g, and 11^ are of the 
form n = n s + n n , g = M(n„v n + n s v s ), and 11^ = P% + M(n n v ni v nj + n s v si v sj ), where n s 
and 77, n are the superfluid and normal fluid number densities, and P is the pressure. Denoting 
the equilibrium quantities by the subscript " " , we have v no = v so = 0, VP D + n o V0 = 0, 
and V(/i + 0) = 0. Using the Gibbs-Duhem relation <i/i = —adT + dP/n, where a = s/n 
is the entropy per particle, these equilibrium conditions imply VT G = 0, and hence 

(i) 

where we have made use of the Maxwell relation {da /8P)t = n~ 2 (dn/dT)p. 

Denoting the deviation of any quantity x from its equilibrium value x Q as Sx = x — x , 
the hydrodynamic equations can be linearized about the equilibrium solution and written 
as I : 5h= —V - (n G v s +n no w), II :n Q v s +n no w = — (5nV0 +V5P)/M, III : v s = -^V{a a 5T 
— — ), IV : <5s = —V ■ (s D v n ), where w = v n — v s . Using eq.([^), I and II imply that 



M8n = V • 



n Q V f — ] — 5Tn V(T 



- A. (2) 



Again using eq.(|]), II and III, we have w = — VffT. By noting that s = cr c ri + ra cr, it 
is easy to show from IV that b = — v n ■ Va — • (n so w), and hence 

Ma = -(Va ) 2 5T + -V • ( ^^V5t) + Va • V ( — )= B. (3) 

n Q \ n no ) \ n Q J 

Expressing all quantities in terms of 5T and 8P, Eqs.(Q) and form a closed set : (|^) T 5P+ 
{jj£) P Sf= A/M, (§p) T 5P + (^) p 5f = B/M. The solutions of these two equations are the 
hydrodynamic modes of the system. To find them, we need to evaluate the thermodynamic 
quantities in these equations. 

Thermodynamics : The thermodynamics of a trapped dilute Bose gas has been worked 
out within the local density approximation (LDA) by Chou, Yang and Yu (CYY) 0. As 



pointed out by CYY, LDA is valid if (a) e = Tiuo/kBT « 1, and (b) A >> a, where 
A = / MksT is the thermal wavelength and a is the s-wave scattering length. These 

conditions are satisfied over a very wide range of temperatures above and below T c for 
large clouds. As discussed in CYY, LDA describes the physics over scales greater than 
a-r = ^Jh/MuT, and all structures on the scale of ar and smaller are shrunk to a point. 
Since the typical width of the interface between the condensate and the normal gas is less 
than ax [0, it is treated as a surface (say, at r = r*) within this scheme. The density is 
continuous at r* but its slope is not ||. Because of this discontinuity, it is necessary to find 
the boundary conditions relating the solutions inside and outside r*. The identification of 
these boundary conditions is the key to our numerical approach, and will be addressed in 
Numerical Methods below. 

Further simplification can be made within the temperature range a\ 2 n <C 1 (denoted 
as condition (c)), where all thermodynamic quantities have been worked out in the classic 
work of Lee and Yang ||. Condition (c) is more restrictive than (a) above. However, it 
still covers a wide range of temperatures. (For a gas of 87 Rb with N ~ 10 6 in a trap with 
lot ~ 10 3 sec _1 , the condition (c) is satisfied over the range 0.6 < T/T c < 1.2 that we 
studied). The coefficients in Eqs.(§) and (|3|) can be calculated in a straightforward manner 
from ref. || and ||. We shall not present the details here for length reasons. Instead, we 
outline the procedure and give the final expressions : 

(i) To determine the density profile at temperature T Q below T c , we first specify the chemical 
potential \x a . This immediately determines the size of the condensate droplet r* through 
the relation \i a = <f>(r*) + 2gn c (T ), where n c (T) = A~ 3 g 3 / 2 (1). (ii) The region r < r* 
consists of both the condensate and the normal components, with n no {r) = n c (T ), and 
n so{ r ) — 5 ,_1 [0( r *) — 0( r )]- The region r > r* consists of only the normal component, with 
n Q (r) determined self consistently from the relation n(r) = X~ 3 g 3 / 2 (C(r)), where ln£(r) = 
P[p, Q — 4>{r) — 2gn (v)]. The quantities iV and fi Q are related though the constraint iV = 
4-7T / drr 2 [n so (r) + n no (r)]. This relation combined with the condition \i = 2gn c (T c ) gives T c 
as a function of N. 
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(iii) When a\ 2 n <C 1, the Hemholtz free energy is /(n(r),T) = —ksTX 3 g 5 / 2 (l) + 
(g/2)[n{v) 2 + 2n{v)n c {T) - n c (T) 2 ] for r < r*, and /(n(r),T) = -/^TA^MCM) + 
fcBTn(r)lnC(r) + gn(r) 2 for r > r* |J. From these two equations, it is straightforward 
to calculate all the thermodynamic quantities needed in Eqs.(^) and (|3|): for example, 
P = —f + n(d f I 'On)? and a = s/n = —n~ l (df/dT) ni and so forth. 

Numerical Method : Since Eq.(^) and (f|) are smooth inside and outside r*, it can 
be solved in each of these regions in a standard way by discretizing it on a grid which is 
made finer as r — > r*. Because some coefficients of these equations are discontinuous at r*, 
(as explained before), the solution of these equations must be understood mathematically 
in terms of the following limiting process. Firstly, we note that the hydrodynamic equations 
I to IV are independent of the specific form of the free energy f(T,n). We can then 
imagine solving these equations for a family of free energies f(T,n;r) parametrized by a 
variable r, which changes from the true free energy f true to the LDA free energy f LDA 
in a smooth manner as r, say, varies from to 1. Such a change, of course, amounts to 
gradually collapsing the actual interface into a very thin region centered at r*. During 
the process of collapse, SP and ST are smooth everywhere. The most rapid changes takes 
place at the boundary of the interface (which we denote as r* ± A), while SP and ST and 
their derivatives are smooth in a close neighborhood (r* — e, r* + e) of r* (e << A) during 
the collapsing process. This implies the following (almost unique) way to implement the 
boundary condition : The interface is modeled by three points r* ± A and r*, where A is 
now the grid spacing. Both SP and ST as well as their derivatives are continuous at r*, 
while the values of SP and ST at r* + A and r* — A are determined by their solutions inside 
and outside r*. Thus, if sharp changes ever occur in the solution, they will only occur at 
r* ± A and nowhere else. 

As we shall see in A below, our numerical method reproduces a set of non-trivial modes 
which can also be obtained analytically. These modes exist in each angular momentum 
sector and are independent of the validity of LDA. The fact that our calculations reproduce 
the exact results for each I shows that the boundary conditions have been implemented 
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correctly. 

Main Results : Because of spherical symmetry, SP and ST can be decomposed as 
(5P(r),8T(r)) = J2e,n,m(.fiPtn(r),STe n (r))Ye m (r), where I is the angular momentum and n is 
the radial quantum number. Each (£, n) mode is 21 + 1 fold degenerate. We have performed 
calculations for 87 Rb (a = 58. 2 A) for the number of trapped atoms N varying from N = 10 4 
to 10 8 . The general features of all these systems are identical. For all cases studied, we 
find that the frequencies of all first and second sound are above the trap frequency ujt H, 
which is to be expected if the system is viewed as a mechanical system with internal degrees 
of freedom in a harmonic potential. For concreteness, we present the results for iV = 10 6 
particles in a trap, with u>t/2tt = 200Hz, over the range 0.6 < T/T c < 1.2 : 

(A) First sound : These modes exist both above and below T c . They are in phase 
pressure and temperature oscillations that extend over the entire cloud, and SP has one 
more node than ST. The frequencies of these modes (denoted as ou^ ) are shown in Fig.l 
for I = 0, 1,2 and n\ = to 4, where n\ counts the number of nodes of SP in the radial 
direction. While {u;^ } change with temperature, their variations are small compared to 
ujt- The eigenfunctions of the [I = l,ni = 2) mode at T = 0.84T CO are shown in Fig. 2a. 
They extend over the entire cloud - a feature common to all first sound modes above and 
below T c . 

The n\ — modes are special. They are isothermal modes of the form SP(r, r) = 
n (r)r ^ m (r), ST = 0, with ujfl = lu t VI. They are also "universal" in the sense that they 
are independent of interaction and statistics. These results emerge from our numerical so- 
lutions but can also be obtained analytically from Eqs.(|2|) and (^) . With ST = 0, using the 
equilibrium relations given by Eq. ([]]), Eqs.(0) and (|3D can be shown to yield V 2 (SP/n ) = 0, 
and df(SP/n ) = — V0- V(SP/n ), which has the solution given above. From the hydrdody- 
namic equations I to IV, it is also straightforward to show that for these isothermal modes, 
v n = v s below T c , and V ■ v„ = both above and below T c . 

The [I = 0, n\ = 1) mode is also special. Above T c , it is a uniform temperature oscillation, 
V5T = 0, but with ST ^ 0. This mode is "non-universal" because it depends on interactions. 
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The interaction effect, however, is sufficiently weak so that ujq~1 is very close to 2ujt above T c . 
It is also straightforward to show that V • v n is constant but nonzero for this mode. These 
results can be established analytically using LDA and also emerge as part of our numerical 



solutions. Below T c , 5T is no longer uniform, and V • v n is not a constant [10|. All the other 
sound modes (£, n\ 7^ 0) are non-isothermal. 

(B) Second Sound : These modes only exist below T c . The frequencies of these modes 
for (i = 0, 1,2) and (n 2 = 0, 1,2) are shown in Fig.l. It should be stressed that the second 
sound frequencies do not merge into the first sounds frequencies asT — > T c . To illustrate this 

clearly, we plot Jf\ 2=x near T c (for £ = to 2) as a function of particle number N in Fig.3. 

(2) 

While uj£n 2 =i changes with N, the first sound frequencies (not shown in Fig.3) typically vary 
by about 2% of uj? in the same range of N. The eigenfunctions of the (i = 1, n 2 = 2) mode 
are shown in Fig.2b. An enlarged structure of the interface of this mode at r* is shown in 
Fig. 2c. The second sound modes have the following common features : 
(1) 8P and ST are "out of phase" inside the condensate and become "in phase" as they leak 
out into the normal region. The leakage reduces to zero as T — > T c . The quantum number 
n 2 counts the number of nodes of 5P or 5T inside the condensate. (2) The wavelengths of 
the oscillations shrink as r — ► r*. This can be understood simply from LDA by recalling that 
the second sound velocity c 2 of a homogenous dilute Bose gas is proportional to ■Jn so . The 
wavelength 2ixk~ 1 is then 2%c 2 /uj oc ^/n so . Since u ~ ujt in our case, and n so (r) vanishes as 
r — > r*, the local wavelength shrinks as r — > r*. (3) The n 2 = modes are different from 
all other n 2 7^ modes. Firstly, except very clost to T c , all ^ 2 2 2= o increase as T decreases, 

(2) 

whereas all <jo\^ 2 , q have opposite behavior, (see Fig.l). Secondly, while 5T and 5P are out of 
phase for all second sound modes (i.e. n 2 = and n 2 ^ 0), the sign of v s ■ v n for a particular 
mode depends on position. In particular, v s ■ v n of the n 2 = modes is actually positive 
(i.e. in phase) almost everywhere inside the condensate instead of negative 0, whereas it 
can be positive or negative for the n 2 ^ modes. (The radial components of v s and v n for 
the n 2 7^ modes are out of phase in most regions in the condensate, so are their tangential 
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components. However, the in-phase and out-of-phase regions of these two components do not 
coincide.) This shows that unlike the second sound modes in homogenous systems, which 
are characterized by either out of phase (5P, ST), or out of phase (v s , v n ) oscillations, the 
correct characterization of the second sound modes in the trap is the out of phase 5P and 
5T oscillations, not the out of phase v s and v n oscillations ||. (4) Near the center of the 
cloud, the ratio £ = \n n v n /n s v s \ is about 0.2 to 0.3 for the modes studied. This is very 
different from 4 He, where the normal current is essentially cancelled by the supercurrent, 
i.e. £ = \n n v n /n s v s \ ~ 1 [§]. That £ is between 0.2 and 0.3 can be understood in terms 
of LDA. From the work of Lee and Yang [Fiji, one finds that £ = ^■(a/X)(gy 2 (l)/g^,/2(l)) 



for the homogenous dilute Bose gas, which is around 0.3 for the temperature range studied 
0. (5) In terms of dimensionless quantities [SP,ST] = [5P/(n (r)k B T ),5T/T ], we find 
that near T c , 5T/5P 3> 1 for all second sound modes, whereas SP ~ 5T for the first sound 
modes. 

Comparison with the JILA data: Examining the JILA data on the sound modes 
of 87 Rb with ~ 2 x 10 3 atoms, we find surprising consistencies with the behaviors of the 
larger systems that we studied : (a) An (m = 0) mode with frequency ~ 2ut was observed 
for all T above T c |2[|. The analog of this mode in a spherical trap is the (£ = 0,rii = 1) 
first sound mode, which also has frequency ~ 2u>t for all T above T c . (b) Below T c , the 
frequency of the (m = 0) mode falls from about 2ut to 1.85^ as T decreases from 0.9T CO to 
0.5T CO p|. The first and second sound analogs of this mode below T c are the (£ = 0, ni = 1) 
and the (£ = 0,n2 = 1) modes respectively. The observed behavior matches well with the 
(£ = 0,n 2 = 1) second sound mode, which drops from about 1.9u>t to 1.5u T as T decreases 
from 0.8T CO to 0.6T CO as seen from Fig.l. (c) An (m = 2) mode was also observed below 
T c 0. Its frequency decreases from 1A5ut to 1.25u;t as T increases from 0.5T CO to 0.85T co . 
The second sound analog of this mode in spherical trap is the (£ = 2, n 2 = 0) mode, which 
also drops from about IAujt to about 1.35co>t as T increases from 0.5T CO to 0.85T co . 

While we do not expect perfect agreement of our results with the JILA observations 
because of the difference in trap symmetry and particle numbers, the qualitative and 
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quantitative consistencies over the temperature and angular momentum range mentioned 
above are striking. The above discussions suggest that the (m = 0) and (m = 2) modes 
observed below T c || are both second sound modes. It is not clear at present why the 
first sound modes do not appear with great prominence below T c . Whether it is due to the 
way that the modes are excited or due to the fact that density oscillations below T c might 
contain a large second sound component because of the large temperature fluctuations in the 
second sound modes (as mentioned in Discussion (5) above) will be studied later. To clearly 
identify the nature of the sound modes, it is necessary to experimentally investigate larger 
number of modes so as to have more consistency checks with the hydrodynamic predictions. 
We hope that this work will stimulate and provide guidance for future experiments. 
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would specially like to thank M.C. Cross for hospitality. This work is supported by NSF 
Grant No. DMR-9705295. 



9 



REFERENCES 

[1] M.H. Anderson, etal, Science 269, 198 (1995). C.C.Bradley, etal, Phys. Rev. Lett. 
75, 1687 (1995). K.B. Davis, etal, Phys. Rev. Lett. 75, 3969 (1995). 

[2] D.S. Jin, etal, Phys. Rev. Lett. 78, 764 (1997). 



[3] Recently, Zaremba, Griffin, and Nikuni ( |cond- mat / 9 70 5 1 34|) have derived the hydrody 



namic equations microscopically in the Popov approximation and have found separate 
conservation of superfluid and normal densities. They have used a variational method 
to obtain the (£ = 1,W2 = 0) second sound mode. They find that the normal and su- 
percurrents are of out phase and cancel each other, with frequencies falling below the 
trap frequency as T — > T c . The characters of the second sound we have found are very 
different from theirs, (see Main Results and Discussions (3) and (4) below). 



[4] CM Kavoulakis, C.J. Pethick, and H. Smith, fcond-mat/9710130|) have shown that 10 7 



atoms are needed to reach the hydrodynamic limit near T c for the Colorado trap. 
[5] Chap. I, Sec. 4 of S.J. Putterman, Superfluid Hydrodynamics, North Holland, 1974. 
[6] T.T. Chou, C.N. Yang and L.H. Yu, Phys. Rev. A 53, 4257 (1996). 
[7] See Fig. 7 of S. Giorgini, L.P. Pitaevskii, S. Stringari, |cond-mat / 97040 14 . 



[8] T.D. Lee and C.N. Yang, Phys. Rev. 112, 1419 (1958). 

[9] For a given grid spacing, all eigenstates except a set denoted as {w)„ 2=0 } below are 
accurate upto one part in 10 4 , while those of the set are accurate upto one part in 10 3 . 

[10] Both w£i 0ni= Q and w^ 0ni=1 modes were also found by A. Griffin, W. Wu, and S. 
Stringari, Phys. Rev. Lett. 78,1838 (1997) for the ideal Bose gas above T c . Interaction 
effects and the behavior of these modes below T c were, however, not investigated. 

[11] T.D. Lee and C.N. Yang, Phys. Rev. 6, 1406 (1959). Also see A. Griffin and E. Zaremba 
( |cond-mat/9707058|) . 



10 



Figure 1: The frequencies ^o\, ni and oj\^ 2 are represented as open and filled symbols resp. 
They are plotted as a function of T/T co , where T c o is the transition temperature of the ideal 
Bose gas in the trap. The dotted line at 0.917T CO indicates the critical temperature T c of 
the interacting model. For the £ ^ modes, r = is not counted as a node. As far as we 
can tell, c^ =0 ,„ 2 =o = 0. 

Figure 2a: The eigenfunctions of ni=2 (marked as A in Fig.l) : To magnify the features 
of the first sound, we have multiplied 5P and 5T by r 2 in Fig. 2a. 

Figure 2b: The eigenfunctions of u;^ =2 (marked as B in Fig.l) : These functions are 
not multiplied by r 2 as those in Fig. 2a because their features are sufficiently clear. 

Fig. 2c : The detailed structure in Fig.2b near r*\ The filled point indicates the location 
of r*, and A is 0.005ay. The sharp change of slope at r* ± A is expected as our boundary 
condition is meant to simulate the collapsing process of LDA at r*. 

Figure 3: The iV dependence of uf^ 2=l modes (£,n 2 = 1) near T c . The temperature is 
chosen so that r* = ax- 
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